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Majorana fermions, originally proposed as elementary particles acting as their own antiparticles, 
can be realized in condensed-matter systems as emergent quasiparticles, a situation often accompa¬ 
nied by topological order. Here we propose a physical system which realizes Landau levels - highly 
degenerate single-particle states usually resulting from an orbital magnetic held acting on charged 
particles - for Majorana fermions. This is achieved in a variant of a quantum spin system due to 
Kitaev which is distorted by triaxial strain. This strained Kitaev model displays a spin-liquid phase 
with charge-neutral Majorana-fermion excitations whose spectrum corresponds to that of Landau 
levels, here arising from a tailored pseudo-magnetic held. We show that measuring the dynamic 
spin susceptibility reveals the Landau-level structure by a remarkable mechanism of probe-induced 
bound-state formation. 


In 1937 E. Majorana suggested [1] that a real wave- 
function can describe spin-1/2 quantum particles that 
are their own antiparticles. Consequently, these particles 
are charge-neutral and decouple from the electromagnetic 
field, making them hard to detect. While no observation 
of Majorana fermions as elementary particles has been 
reported to date, their possible realization as emergent 
quasiparticles in condensed-matter systems has attracted 
enormous attention. Specifically, it was realized that su¬ 
perconductors provide a natural habitat, where Bogoli- 
ubov quasiparticles at zero energy have properties of Ma¬ 
jorana fermions [2]. Moreover, a seminal work of Ki¬ 
taev [3] demonstrated that dispersive Majorana fermions 
can emerge as effective degrees of freedom in spin-liquid 
phases of frustrated quantum magnets [4]. This has trig¬ 
gered an intense search for materials which come close to 
realizing Kitaev’s spin-liquid model on the honeycomb 
lattice, with A 2 lr 03 (A=Na,Li) (Refs. 5-9) and a-RuCla 
(Ref. 10) currently being the best candidates. 

Here we show that the properties of emergent Majo¬ 
rana fermions in quantum magnets can be engineered 
by controlled lattice distortions. We use spatially inho¬ 
mogeneous exchange couplings, generated by applying a 
suitable strain pattern, to transform the linear Dirac-like 
dispersion of low-energy Majorana fermions in the Ki¬ 
taev model into a sequence of pseudo-Landau levels. The 
idea of strain-induced Landau levels was in fact first the¬ 
oretically proposed [11] and then subsequently realized 
[12] for electrons in two-dimensional graphene. Our work 
adapts this concept to the charge-neutral fractionalized 
quasiparticles of a topological spin liquid. We compute 
the dynamic spin susceptibility which turns out to dis¬ 
play sharp excitations at isolated energies reflecting the 
Landau-level structure of the strained spin liquid. 

Model. The Kitaev model [3] describes quantum spins 
1/2 on a honeycomb lattice subject to spin-anisotropic 
Ising (or “compass”) interactions. The Hamiltonian, gen¬ 


eralized to spatially varying couplings, reads 

Wk = - E - E ■’tNN - E ■'.A'"} (1) 

(b>x {ij)y {ij)z 

where (t“ are Pauli matrices, and {ij)a denotes an a bond 
as shown in Fig. 1, with a = x,y, z. This model displays 
an infinite set of constants of motion, which can be in¬ 
terpreted as Z 2 fluxes through the closed loops of the 
lattice. Upon representing each spin a by four Majorana 
fermions b^, IP, and c, with af = ibfci, the model (1) 
takes the form 

"Hu = iX JijUijhCj, (2) 

iu) 

where Uij = i6“‘^6“‘^ and Uij = —Uji. The operators 
Uij represent conserved quantities with eigenvalues of 
Uij = ±1 which determine the values of the Z 2 fluxes. 
Hence the Hamiltonian Hf, (2) represents a nearest- 
neighbor hopping problem for the c (or “matter”) Majo¬ 
rana fermions which are coupled to a static Z 2 gauge field. 
Its exact solution can be written in terms of canonical- 
fermion excitation modes with non-negative energies Cm- 
For homogeneous and isotropic couplings, J/)- = J, it de¬ 
scribes a Z 2 spin liquid with matter-fermion excitations 
displaying a gapless Dirac spectrum. 

For our case of inhomogeneous couplings we shall solve 
the Kitaev model numerically on finite-size lattices. To 
this end, Eq. (2) is re-written into the following bilinear 
Hamiltonian for the matter Majorana fermions 

^ (c5 cl) ( 3 ) 

where M is a.n N x N matrix with elements = J^jUij 
and Ca(b) is a vector of length N of Majorana opera¬ 
tors on the A{B) sublattice, and N the number of unit 
cells. We construct the matrix M for a given set of cou¬ 
plings {J^j} diagonalize the problem via a singular- 
value decomposition which yields the excitation energies 
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FIG. 1: Lattice setup: Honeycomb-lattice flake with N = 675 
unit cells where bonds with different color correspond to Ising 
interactions J“ of different spin components a. The distorted 
lattice visualizes the displacement according to Eq. (4) arising 
from triaxial strain; the unstrained hexagonal flake is shown 
in the back in grey. Longer (shorter) bonds correspond to 
weaker (stronger) exchange couplings Jj“ . 


Cm and the eigenmodes. From these we calculate the lo¬ 
cal density of states (LDOS) as the local spectral density 
of the c fermions, for details see supplement [14]. 

Majorana Landau levels. In the context of graphene, 
it has been theoretically shown [11, 15] that a spatial 
modulation of hopping energies mimics the effect of a 
vector potential in Dirac fermion systems. Near the 
Dirac energy, this emergent vector potential can be ex¬ 
pressed through the strain tensor Uij as A (x i(uxx — 
Uyy, — 2uxy)^, with opposite sign for electrons belong¬ 
ing to the two Dirac cones (or “valleys”) located at 
K = (47r/(3-\/3), 0) and K' = (—47r/(3-\/3), 0) (with the 
lattice constant oq set to unity), such that time-reversal 
symmetry is preserved. If the resulting pseudo-magnetic 
field, B = rotH, is sufficiently homogeneous - applying 
e.g. to triaxial strain patterns as shown in Fig. 1 - it can 
induce single-particle pseudo-Landau levels very similar 
to Landau levels in a physical magnetic field. The system 
then displays a so-called quantum valley Hall effect [11], 
Le., it combines a chiral and an antichiral quantum Hall 
effect at the valleys K and K', respectively. For triaxial 
strain the displacement vector is given by [11, 13] 

U{x,y) = C{2xy,x^-y^) (4) 

where C (measured in units of I/oq) parameterizes the 
distortion, and Uij = {diUj-\-djUi)/2. This displacement 
yields - to first order in C - a homogeneous pseudo- 
magnetic field whose strength is proportional to C. 

Adapting the idea of strain-induced artificial gauge 
fields to spin systems, we study the Kitaev model (1) with 
spatially modulated exchange couplings. Given that the 
Majorana-fermion hopping matrix elements in Eq. (2) are 
given by the exchange couplings J“-, we choose 


j“ = r 


1 - ^(l^ijl/ao - 1) 


(5) 


where we calculate the distance Sij = Ri + Ui — Rj — Uj 
using the U{x,y) from Eq. (4) evaluated at the lattice 
positions Ri of the undistorted honeycomb lattice. The 
factor fd encodes the strength of magnetoelastic coupling, 
and the parameter (dubbed “strain” below) C = (7/3 will 
enter our simulations as a measure of the modulations of 
the J“-. Throughout the paper we choose /3 = 1; larger 
values of /3 (hence smaller (7) reduce non-linearities in 
the strain pattern but are less realistic for transition- 
metal oxides. We note that, for a real material, Eq. (5) 
represents a linear approximation to the full dependence 
of the exchange constant on the bond length [13]. 

The inhomogeneous Kitaev model (1) with couplings 
given by Eq. (5) is expected to display a Z 2 spin-liquid 
phase with an extremely unusual excitation spectrum: 
While the Z 2 fluxes remain static, the matter Majo¬ 
rana fermions will form highly degenerate Landau levels. 
Given that a physical magnetic applied applied to a Ki¬ 
taev magnet has an entirely different effect - it induces 
flux dynamics and leads to a gapped matter spectrum [3] 
- strain represents a unique way to generate Majorana 
Landau levels. 

This is well borne out by our numerical results, ob¬ 
tained for finite-size Kitaev systems of hexagonal shape, 
Fig. 1, with open zigzag boundaries and lattice sizes up 
to 2N = 15 000 sites. We focus on the case of isotropic 
couplings and choose = J as energy unit. We limit 
the strain C to be smaller than (7inaxj the latter being 
defined as the largest C where all > 0 [14]. For 
(7 < (7max we find that the ground state is located in the 
flux-free sector where all Uij can be chosen to be -1-1. We 
monitor the LDOS of the matter fermions [14] which is 
shown in the top panels of Fig. 2: In the unstrained case 
(7 = 0 this is the familiar honeycomb-lattice DOS [panel 
(a)], with p{ui) cx w at low energies, with the difference 
to graphene that only half of the spectrum correspond¬ 
ing to w > 0 is realized due to the Majorana nature of 
the matter fermions. Importantly, the results for finite 
strain show clear Landau-level peaks at low energies high¬ 
lighted by the black arrows in panels (d) and (g). In par¬ 
ticular, the Landau-level energies display the scaling 
£n oc Vu(7 with the lowest Landau level (LLL) located 
at zero energy, Fig. 3, characteristic of honeycomb-lattice 
Dirac fermions [16]. 

Dynamic structure factor. In order to detect the 
Landau-level structure, we propose to measure dynamic 
spin correlations 

( 6 ) 

whose Fourier transform is - in a magnet - accessible by 
neutron scattering experiments. We restrict our atten¬ 
tion to zero temperature where {uj) is proportional 
to the imaginary part of the dynamic spin susceptibility. 
The evaluation of the dynamic spin correlations in the 
Kitaev model has been first discussed in Ref. 17. Specifi¬ 
cally, the application of the operator df = ibfci changes 
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FIG. 2: Simulation results: Histogram of the matter-fermion LDOS in the flux-free (top, blue) and two-flux sectors (middle, 
green), together with the dynamical spin correlator S^ij){Lo) with S(^ij) = Su -f Sjj + 2Sij (bottom, purple), for the cases (a-c): 
unstrained (C = 0), (d-f): moderate strain (C = 1 • 10“®), and (g-i): large strain ((7 = 4- 10“®). All quantities are measured 
in the center of a sample with 2N = 15 000 sites. The S{uj) plots employ a Lorentzian broadening of width 7 ; moreover the 
energy axis has been shifted with respect to the LDOS plots by the amount of the local flux gap[14] Aij (yellow region). Panel 
(c) also shows S{uj) for the infinite homogeneous system [17]. While the LDOS of the flux-free sector [panel (d) and (g)] nicely 
shows the emergence of Majorana Landau levels (black arrows), the two-flux sector features bound states in the gaps between 
the Landau levels [panels (e) and (h)] - these are detected in as sharp peaks (dashed lines are guide-to-the-eye). 


the Z 2 fluxes in the plaquettes adjacent to the a bond 
emanating from site i. Starting from the ground state 
jo) which is flux-free, the intermediate state |A) is then 
a state with a locally excited flux pair, i.e., the hopping 
problem of matter fermions in the intermediate state in¬ 
volves a local “flux impurity” V = —2i The fact 

that fluxes are static implies that all correlators beyond 
nearest-neighbor sites are strictly zero. 

The numerical calculation of a zero-temperature dy¬ 
namic structure factor Sij{uj) involves two diagonaliza- 
tions, one in the flux-free sector with all u = 1 and one 
in the excited-state flux sector with all rt = 1 except 
Uij = — 1 on the bond involving the measured site(s). 
The eigenvectors in both flux sectors are then used to 
construct the matrix elements entering Sij{(jj). We em¬ 
ploy the single-mode approximation [17], with details 
given in the supplement [14]. 

For the structure factor of the strained Kitaev model, 
a central observation is that the impurity V induces a se¬ 
quence of single-particle bound states which energetically 
he between the Landau levels. This can be nicely seen 
in panels (e) and (h) of Fig. 2 which show the LDOS 
in the two-flux sector, measured on the flipped bond. 
This LDOS displays essentially the same Landau levels 
as the corresponding flux-free system, but shows addi¬ 
tional sharp peaks of even higher weight at energies En 


in between the Landau levels - these peaks arise from 
isolated states localized near the flipped bond. We note 
that related bound states have been reported recently 
in a study of spin excitations in the field-induced non- 
Abelian phase of the Kitaev model [18]. More broadly, we 
recall that in-gap bound states occur frequently in one- 
and two-dimensional systems because the real part of the 
single-particle Green’s function diverges at gap edges. In 
fact, recent theory work on topological band insulators 
has suggested that the generic occurrence of impurity- 
induced bound states in a two-dimensional gapped sys¬ 
tem is a signature of its topological nature [19]. For the 
present case of bond impurities in the triaxially strained 
honeycomb lattice, we have checked that bound states 
between the Landau levels occur indeed generically, i.e., 
independent of the impurity strength [14], hinting at the 
topological character of the matter-fermion sector. 

Remarkably, these Majarona-fermion bound states en¬ 
able a unique detection of the Landau-level structure: 
The largest matrix elements for the correlator (w) 
arise from these bound states in the two-flux sector, such 
that (ui) displays sharp peaks at E^ + A^ where A^- 
is the local flux gap [14, 20]. As can be seen in the lower 
panels of Fig. 2, these bound-state peaks dominate over 
the peaks arising from the Landau levels themselves. 

In Fig. 3 we show the energetic positions of the bound- 
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state peaks in with the local flux gap A^- sub¬ 

tracted, and compare them to the energies of the Landau 
levels found in the matter-fermion LDOS of the flux-free 
system. The plot shows that the bound-state energies 
track the energy dependence e„ oc '/nC of the Landau 
levels, approaching the lower edge of the respective gap 
with increasing n. Given that sharp bound states can 
only arise in gapped systems, the observation of a series 
peaks in with a ^/n energy dependence (after 

subtracting a constant corresponding to Aij) represents 
a direct signature of the peculiar excitation spectrum - 
Landau levels separated by gaps - of the strained spin 
liquid. 

Discussion. Our results raise a number of further is¬ 
sues. First, we note that the simultaneous application of 
strain and of a physical magnetic field has non-trivial ef¬ 
fects, which can be analyzed perturbatively by projecting 
onto the flux-free sector [3]. We have found that the LLL 
of the strained Kitaev model is shifted to a finite energy, 
while otherwise the spectrum of the matter Majorana 
fermions remains qualitatively unchanged. The excita¬ 
tions of the now fully gapped system obey non-Abelian 
statistics similar to the elementary excitations of Kitaev’s 
B-phase [3] subject to a magnetic field. Second, and per¬ 
haps most interesting, is the fate of the Landau-level spin 
liquid when interactions beyond the Kitaev model are in¬ 
cluded [21]. While this question is beyond the scope of 
this work, it is likely that the large degeneracy of the 
many-body ground state will be lifted in favor of po¬ 
tentially exotic states driven by interactions between the 
Majorana fermions, akin to the physics of the fractional 
quantum Hall effect. 



FIG. 3: Energy positions of the Landau levels as observed 
in the matter-fermion LDOS of the flux-free sector, shown to¬ 
gether with the bound-state energies En in the two-flux sec¬ 
tor; the En are identical to the peak positions in the correlator 
with the local flux gap Aij subtracted [except for the 
lowest peak which originates from the LLL; hence our num¬ 
bering for En starts with n = 1]. The data correspond to 
that of Fig. 2(g,i) with C = 0.004. The lines are fits of £n 
to a ^/n dependence; both £„ and Sn+i are shown for visual¬ 
ization purposes such that the Landau-level gap is visible for 
each n. 


Realizations of our system appear possible using 
strained thin films of honeycomb-lattice iridates, where it 
has been proposed that homogeneous pressure or strain 
may drive the material into a Kitaev spin-liquid regime. 
Landau levels can then be generated using shapes result¬ 
ing from triaxial strain, circular arcs [22], or nanobubbles 
[12]. It has to be kept in mind, however, that inhomoge¬ 
neous strain may induce additional interactions as it de¬ 
forms the nearly 90-degree Ir-O-Ir bonds. A second possi¬ 
bility is offered by advances in designing artificial molecu¬ 
lar structures on surfaces, where the electronic structure 
of strained graphene has already been demonstrated [23]. 
A third option is by simulating the Hamiltonian (1) us¬ 
ing ultracold atomic gases. Phase plates [24] or digital 
mirror devices [25] can be used to directly project the 
distorted lattice onto a 2D cloud of atoms, and proposals 
for realizing artificial spin-orbit coupling have been made 
[26]. 

In summary, we have demonstrated how to engineer 
a novel state of matter, where the emergent degrees of 
freedom of a fractionalized spin liquid - charge-neutral 
Majorana fermions - display a Landau-level structure of 
excitation energies. We have shown that this Landau- 
level structure can be efficiently detected in measure¬ 
ments of dynamic spin correlations, thanks to a mech¬ 
anism of probe-induced bound-state formation. 

We thank J. Knolle, Y. M. Lu, T. Meng, M. Neek- 
Amal, F. Peeters, and F. Zschocke for discussions. This 
research was supported by the DFG through SFB 1143 
and GRK 1621 as well as by the Helmholtz association 
through VI-521. 
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I. EXACT SOLUTION OF THE KITAEV 
MODEL 


The spin-liquid model of Kitaev^ describes spins 1/2 on 
a honeycomb lattice which interact via bond-dependent 
nearest-neighbor Ising spin exchange, 

= - E E E (SI) 

(b>x {ij)y {ipz 

with af being Pauli matrices. For the homogeneous 
isotropic case, J“- = J, the Kitaev model possesses a 
discrete Z 3 symmetry of a combination of 120 degree 
real-space and simultaneous spin rotations of the form 
(jv ^ cry —>• ( 7 ^, and ^ . As usual, the two 

sublattices of the honeycomb lattice will be labelled A 
and B. 

Following Kitaev^ one rewrites the spin Hamiltonian 
(SI) in terms of four Majorana fermions per site, 6 ^, 

6 ^, and c, with uf = ibfci. In the new operators, the 
Hamiltonian (SI) becomes 


"H* = i E "^5 “b c* Cj (S2) 

(ii)a, 


with the “bond operator” Uij = ibfb'j which is antisym¬ 
metric under exchange of i and j (and i taken from the 
A sublattice). The problem can be solved by noting that 
all Uij are constants of motion: they commute with the 
Hamiltonian "H* and with each other, their eigenvalues 
being Uij = ±1. Every particular choice of a given set 
{uij} reduces the Hamiltonian (S2) into a billinear which 
is readily solved, 


n 


u — 



0 M 
0 



(S3) 


which can be combined into complex fermions 


— 2 (^A,m + i c'b,™) 
Using the matrix identity 


(S5) 


/ 0 USV'^\ /[/ 0 W 0 5*/ 

\-VSU^ 0 

we can rewrite the Hamiltonian (S3), 


0 

[ 0 

(S 6 ) 










0 

-s 




N N 

= i ^ ^ ' Cm ~ 1 ) ■ (^'^) 

m—1 m—1 


The matter-fermion excitation energies are thus given by 
the singular values Cm of the matrix M ; the ground-state 
energy for a given choice of {uij} is Ao,„ = -J2m^rn- 
Following Kitaev we will refer from now on to the choices 
{uij} as flux sectors. The case where all Uij = +1 is 
referred to as the flux-free sector. The case where all 
Uij = - 1-1 except a single = — 1 is referred to as 
a single-flux sector if the bond (fj-v) is located at the 
boundary of the sample and as two-flux sector otherwise. 

To monitor the eigenmodes of the matter Majorana 
fermions, we can calculate their (global) density of states 
(DOS), p{u)) = J2n ~ 2en)- To resolve spatial struc¬ 
ture, we also consider the corresponding local density of 
states (LDOS) which is given by 


For a honeycomb lattice with 2N lattice sites, M is a real 
NxN matrix with elements Mj^. = Jp,Ujk and the vector 
ca(b) contains N matter Majorana operators of the A 
(B) sublattice. Using a singular value decomposition one 
finds M = USV"’" with orthogonal NxN matrices U 
and V and the positive semi-definite diagonal matrix S 
with elements e*, i = 1,... ,N. Denoting the Majorana 
fermions on sublattice A (B) in unit cell r as CA(B),ri the 
Majorana eigenmodes are given by 

^A,m ~ ^ ' ^m,r^A,r and Cg „^ = ^ ' Vjy, j.CB,r (S4) 


Pj{uj) = lj2^njSiu;-2e^) (S 8 ) 

n 

for a site j on the A sublattice; for a site j on the B 
sublattice, Unj is replaced by Vnj where Unj (Uij) are 
the matrix elements oiU (V). 

In our numerical calculations we consider hexagonal 
flakes with open boundaries (see Fig. SI), consisting of r 
rings, with the number of unit cells being N = 3r^. Most 
runs are performed for r = 50, yielding 2N = 15 000 
spins. 
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FIG. SI: Illustration of a hexagonal shaped honeycomb flake 
consisting of r “rings”. Shown are r = 1 (white), r = 2 
(yellow + white), r = 3 (red + yellow + white). The number 
of unit cells is = 3r^ corresponding to 2N spins. 

II. DYNAMICAL SPIN CORRELATIONS 

The dynamical spin correlation function, Eq. (5) of 
the main paper, can be expressed in terms of Majorana 
fermions and yield the following Lehmann representation 
at zero temperature^’^, here written separately for the 
on-site correlator 

5“^(a;) = 2^^ I (Mol c, |A) fd{cu- [Ex-Eo))d^p (S9) 

A 

and the off-site correlators 

S^fico) = ^ (Mol c, |A) (A| c, |Mo) (SIO) 

A 

X (5(cu {E/_\ E/o))^Q/3^(ij)a 5 

where is non-zero only if i and j are nearest neigh¬ 

bors connected by an a bond, i.e., Sij vanishes beyond 
nearest neighbors. In both cases |Mo) is the matter- 
fermion ground state in the ground-state flux sector, 
assumed to be flux-free, and runs over all matter- 
fermion states in the two-flux sector. Eq and E\ are the 
corresponding many-body energies. The prefactor = 
{ —1, i, — i} depending on the spin component. Note that 
we can safely ignore the fermion parity condition^’^ due 
to open boundaries, such that |Mo) has no excited matter 
fermions. For the intermediate state |A) we will restrict 
ourselves to states with a single matter fermion (single¬ 
mode approximation); these states contribute more than 
97% of the spectral weight in the homogeneous case^. 
With strain applied, we can expect this approximation 
to be even better, because the main weight is contained 
in the bound-state peaks whose contributions are cap¬ 
tured exactly in the single-mode approximation. 

In the homogeneous isotropic Kitaev model, the 
momentum-dependent dynamic spin structure factor 

- accessible by 

neutron scattering - displays a continuum of spectral 
weight above the energy of the flux gap, A = 0.26 J, with 
little momentum dependence^. This reflects the fact that 
a spin flip decays into a matter and a flux excitation, the 
latter having a momentum-independent energy. 


For our inhomogeneous Kitaev model, the computa¬ 
tion of the spin correlator proceeds in the following steps: 

1. We construct the matrix M in (S6) for the flux-free 
sector (i.e., all Uij = -|-1). This matrix is real and 
of dimension N x N, its entries require to specify 
the values of J“ and the strain C. 

2. We perform a singular value decomposition (SVD) 
into M = USV"’" where U and V contain the 
normalized singular vectors and are orthogonal, 

= U~^ and = V~^. The diagonal matrix S 
contains the singular values - these are identical to 
the non-negative part of the spectrum of an equiv¬ 
alent hopping problem of canonical fermions with 
hopping energies Mjk- 

3. We obtain the M' matrix of a particular two-flux 
sector, defined by a single sign-flipped bond (ij) 
(equivalent to multiplying a single matrix element 
in M by —1). Now we also perform the SVD for 
the matrix M' and obtain corresponding U' and 
V matrices as well as the singular values . We 
denote the matter-fermion ground state in this two- 
flux sector by |Ao) and the canonical fermions de¬ 
scribing the corresponding excitations [analogous 
to those in Eq. (S5)] by bm- 

4. We define the matrix 

X = ^(u'^u+ v'^v'^ (Sll) 

and compute its determinant and its inverse, det X 
and X~^, respectively. 

5. For an intermediate state with a single matter- 
fermion excitation, |Am) = |Ao)j matrix ele¬ 
ments for c operators on the A and B sublattices 
are obtained as 

{Mo\cA,^\X^) = V|detX|(C/A-i)^_^ , (S12) 

{Xm\cB,j\Mo) = iV|detA|(FA-i)^._^ . (S13) 

6. Eventually we calculate the correlators (S9) and 
(SIO), with the ground state energy of the flux-free 
sector Eq and with the one-particle energies of the 
considered two-flux sector, with 

Exm = E^q^ + 2e(^) =-Y,et]+ 2e^^ • (S14) 

m' — l 

III. FINITE-SIZE AND BOUNDARY EFFECTS 

The displacement pattern as described by Eq. (3) is 
incompatible with periodic boundary conditions, such 
that open boundaries need to be employed. Then, the 
pseudo-magnetic field in a finite-size system as in Fig. 1 
is not perfectly homogeneous. In addition, low-energy 
edge states occur. 
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A. Maximal strain 

The assumed linear dependence of the exchange cou¬ 
plings on bond lengths, Eq. (4), implies that there ex¬ 
ists a maximal value of the strain, Cmax, such that for 
C > Cmax exchange couplings of the boundary of the 
sample become negative. The value of Cmax decreases 
with system size; concrete values are Cmax = 0.021 for 
2N = 1350 and Cmax = 0.0062 for 2N = 15000. All re¬ 
sults shown in this paper are for strain values C < Cmax- 
We have numerically checked that the ground state of 
the Kitaev model remains in the flux-free sector provided 
that C ^ Cmax- 


B. Landau levels and field inhomogeneities 

As a result of the field inhomogeneities, the clearest 
Landau-level signatures are obtained in the LDOS in the 
center of the system. This has been discussed in detail 
for graphene^, and our LDOS results including finite-size 
effects are similar to that in Ref. 5. 

We note that, despite the inhomogeneities, the pseudo- 
Landau levels display a large approximate degeneracy, 
i.e., the energy of the Landau-level peaks is essentially 
constant as function of spatial position in an extended 
region around the sample center. To illustrate this, we 
show in Fig. S2 the LDOS of Fig. 1 g (main paper), mea¬ 
sured at the sample center, and two corresponding LDOS 
further away from the center. In all three cases the Lan¬ 
dau levels appear at the same energy. Close to the sam¬ 
ple edge the energy position of the Landau levels slightly 
deviates (with the exception of the LLL). We note that 
larger values of j3 reduce non-linearities in the strain pat¬ 
tern which are ultimately responsible for the deviations 
between Figs. S2(a), (b), and (c). 


C. Edge states 

While our system displays edge states at its zigzag 
edges, these are not topologically protected and do not 
energetically connect different Landau levels - the latter 
can be clearly seen in our LDOS results for large values 
of strain C (not shown). This is consistent with the fact 
that the strained system remains time-reversal invariant 
and hence is characterized by a vanishing Chern number. 


D. Flux gap for inhomogeneous couplings 

For the homogeneous isotropic Kitaev model the flux 
gap, be., the energy difference between the matter- 
fermion ground states in the two-flux and zero-flux sec¬ 
tors, is given by A = 0.26J in the thermodynamic limit^. 
Note that the lowest-energy two-flux states are those with 
the two fluxes being located in adjacent plaquettes. 
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FIG. S2: LDOS in the center and off the center for 2N = 
15 000 spins and C = 0.004. (a) is identical to Fig. 2(g), 

measured in the center, (b) and (c) show LDOS measured 10 
and 20 sites away from center, respectively. Red dashed lines 
are guide-to-the-eye for comparing the position of the LLs. 


Importantly, for a finite-size system with open bound¬ 
aries the flux gap depends on the spatial position of the 
flux pair. This requires to introduce a local flux gap A^ 
where ij refers to the bond with = — 1 adjacent to 
the fluxes. It is this local gap A^- which represents the 
lower energy bound for the intermediate state |A) in the 
spin correlator 

In Fig. S3 we show A^ for different values of the strain 
C. Clearly, small and moderate values of C lead to an 
essentially homogeneous gap in the interior of the system 
and a reduced gap near the edges only. In contrast, for 
large C close to Cmax the gap profile becomes rather inho¬ 
mogeneous, with a pronounced maximum in the middle. 
This can be easily understood: The strongly elongated 
bonds near the edges correspond to a reduced exchange 
coupling, and the energy cost of flipping a bond in such 
a plaquette is small. 































































































































































































































4 



FIG. S3: Local flux gap Aij / J. The four panels correspond to 
different values of C: (a) C = 0, (b) G = 0.008, (c) C = 0.016, 
and (d) C = 0.02. The sample contains 2N = 1350 lattice 
sites. 


E. Local vs. zero-momentum spin correlations 

Given the inhomogeneous character of the system, with 
both the pseudo-magnetic field and the flux gap being 
approximately homogeneous in the interior only, one has 
to distinguish measurements in position and momentum 
space. Momentum-space spin correlations are probed by 
the dynamic structure factor S'(g, u). For real-space mea¬ 
surements it is convenient to introduce a bond-local re¬ 
sponse function, 

{u:) = (^) + + 2^5“ (^) (S15) 

with a = aij. The structure factor is related 
to these local correlators by — 0,a;) = 

^/^J2{ij) where the sum {ij) is over all lattice 

bonds. 

The plots of Fig. 2 in the main paper display local ob¬ 
servables in the center of the system where the Landau- 
level signatures are clearest - here isolated bound-state 
peaks are seen in Importantly, the condition of 

locality for the measurement can be relaxed: In Fig. S4 
we show approximants to (1/3) = 0,a;), de¬ 

fined as the sum of all correlators inside a circle of radius 
R in the interior of the sample: 

ri,rj<R 

^ S(ij)iu}) (S16) 

{ij) 

where N' the number of bonds inside this circle. For 
small i?. Fig. S4(a), the result for Sr resembles the single¬ 
bond result. Increasing R leads to a smearing of the 


bound-state peaks - this smearing primarily originates 
from the inhomogeneities in the flux gap. Importantly, 
the bound-state peaks characteristic of the Landau-level 
structure are still clearly visible for R values where the 
circle covers a sizeable fraction of the sample, 20% in 
Fig. S4(c). This shows that neutron-scattering measure¬ 
ments can detect the peak structure by measuring S{q, w) 
provided that the neutron beam is smaller than the sam¬ 
ple size, such that the outer parts of the sample are not 
included in the probe. 


F. Thermodynamic limit 

The displacement pattern used to generate the pseudo- 
magnetic field, Eq. (3), cannot be extended to an in¬ 
finitely large system, simply because bond lengths be¬ 
come negative. In other words, the thermodynamic limit 
of infinite linear size, L/oo —> oo, cannot be taken at fixed 
C. We also note that, with increasing system size, the lin¬ 
earization in Eq. (4) becomes invalid, which is eventually 
signalled by negative exchange couplings at the bound¬ 
ary which occur for C > Cmax- Together, this raises the 
question whether there exists a well-defined limit where 
finite-size effects can be made small. 

An interesting route is to take the combined limit 
L/ao —>■ oo and C —0, keeping CL/gq (as well as the 
magnetoelastic factor /3) fixed. In this limit, the pseudo- 
magnetic field scales to zero as L/qq —> oo, and with it 
the energetic spacing of the pseudo-Landau levels. As the 
ratio of system size and magnetic length is kept fixed, we 
speculate that the system approaches a scaling limit with 
well-defined Landau levels - this will be investigated in 
future work. 


G. Spin correlations at finite temperature 

For an orbital magnetic field applied to graphene, the 
lowest Landau level (LLL) is located at exactly zero en¬ 
ergy. For the present case of the strained Kitaev model, 
the LLL being at zero energy implies a highly degener¬ 
ate many-body ground state. Notably, this degeneracy is 
lifted by finite-size effects - all our mode energies Cm are 
positive. 

For practical reasons, we evaluate the expectation 
value (...) in the spin correlator, in Eq. (5), with the 
unique finite-system ground state \Mq), corresponding to 
strictly T = 0. Importantly, our conclusions apply un¬ 
changed if is calculated instead using a thermal aver¬ 
age over the low-energy many-body states, i.e., at small 
finite temperature T <C ei where ei is the energy of the 
first finite-energy Landau level. The argument proceeds 
as follows: The relevant many-body initial states \M) in¬ 
volve any number of matter Majorana excitations from 
the LLL - this precludes a full numerical evaluation of 
5'“^ at finite T. Importantly, for large systems, sizeable 



























5 





FIG. S4: Correlators Sr{uj) (S16) - representing a finite-area approximants to S{q — 0,oj) - for (a) i?/i?max = 0.1, featuring 
essentially the same peak structure as S^ij){uj) in Fig. 2(i), (b) i?/7?max = 0.2, with visible broadening of the bound-state peak 
due to inhomogeneities, and (c) i?/i?max = 0.4. The system parameters are 2N — 15 000, C = 0.004, and 7?niax is the radius of 
the smallest circle enclosing the entire flake (in unstrained coordinates). The Lorentzian broadening is 7 /J = 0.005. 
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FIG. S5: Wavefunctions of the first two in-gap bound states 
in the two-flux sector for a system with 2N — 5400 spins and 
strain C = 0.0078 where the flux pair has been created in the 
center, (a) First in-gap state with energy En^i « 0.44 J. (b) 
Second in-gap state at En =2 ~ 0.86 J. 


matrix elements (A| Cj \M) are only obtained if the inter¬ 
mediate state contains, in addition to the excitation ef¬ 
fectively created by Cj , the same excitations as the initial 
state \M). This has been demonstrated in Ref. 3 where 
the effect of the ground-state fermion-parity condition^’^ 
on the spin structure factor has been examined: It was 
found that the two-particle continuum, obtained by start¬ 
ing from an initial state with one excited fermion, is vir¬ 
tually indistinguishable from the one-particle continuum, 
obtained by starting from an excitation-free initial state. 
More precisely, the differences between the two results 
for 5'“-^ scale with the inverse system size, the reason 
being that the wavefunctions of the additional excita¬ 
tions are spatially extended. Generalized to our system, 
this implies that all thermal contributions to the low- 
temperature spin correlator will have essentially the same 
peak structure, arising from the bound states between the 
Majorana Landau levels. While small finite temperatures 
will cause both small broadening and a small weight re¬ 
distribution, the peak structure in S(l>j) shown in Fig. 2 
can be expected to be robust. 


IV. PSEUDO-LANDAU LEVELS AND BOUND 
STATES 


As emphasized in the paper, sufficiently large strain 
induces pseudo-Landau levels in the matter Majorana 
sector of the Kitaev model. At the same time, the dy¬ 
namic spin correlator S(^ij){uj) shows sharp S peaks. In¬ 
terestingly, these peaks (with the exception of the low¬ 
est one which can be attributed to the LLL in the two- 
flux sector) do not directly arise from the Landau levels. 
Instead the measurement dynamically generates an im¬ 
purity leading to one localized bound state, causing a 
high-intensity peak in in each of the Landau- 

level gaps. Given that higher Landau levels tend to be 
smeared (because the analogy between strain and mag¬ 
netic field is restricted to the low-energy limit), the spin 
correlator shows a much sharper peak structure than the 
matter-fermion DOS itself. 


A. Bound-state wavefunctions 

In Fig. S5 we illustrate the nature of the probe-induced 
bound states in the two-flux sector, by plotting the prob¬ 
ability distribution of eigenmode m for A and B sites 
given by 

= 1^4/2. (S17) 

The bound-state eigenmodes are highly localized near the 
flipped bond (which is taken to be in the center of the 
sample), which we have checked by computing the in¬ 
verse participation ratio for these modes. The localiza¬ 
tion length increasing with increasing n, which is a simple 
result of the decreasing gap magnitude. 


B. Bound states and topology 

A recent paper® has proposed an interesting connection 
between topological states and local impurity physics, 
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FIG. S6: Energy positions En of the bound states caused 
by a bond impurity as a function of the impurity strength 
Vimp. Note that Fimp = 0 is the clean case and Vimp = 2Jij 
corresponds to a flipped bond (“two-flux sector”). The data 
corresponds to N = 7500 and C = 0.004; the grey horizontal 
lines represent the energies of the pseudo-Landau levels. 

namely that the impurity-induced local in-gap Green’s 
function can be used as diagnostic tool for topology. For 
time-reversal-invariant Z 2 topological insulators in two 
space dimensions a local impurity always causes an in¬ 


gap bound state. This is because the real part of the 
local Green’s function diverges at the lower edge of the 
upper band with opposite sign compared to the upper 
edge of the lower band. This guarantees the presence of 
an in-gap bound state. It was further shown that a topo¬ 
logically trivial insulator does not generally possess such 
in-gap bound states (but may do so in a finite interval of 
impurity strengths). 

As our analysis has revealed the creation of localized 
in-gap states due to bond impurities, the latter repre¬ 
senting a pair of Z 2 fluxes, we adapt the idea of Ref. 6 
to our setting. We consider a bond impurity of general 
strength, V = — iFimpCiCj, in a honeycomb-lattice Ma- 
jorana hopping problem subject to triaxial strain. While 
the impurity value Ijmp = — 2 (measured in units of Jy ) 
corresponds to a flipped bond, i.e., a two-flux sector, here 
we consider arbitrary values of Ijmp G [~15,15], with 
the clean case corresponding to Ijmp = 0. As shown in 
Fig. S6, we find that bond impurities of any strength 
induce bound states between the Landau levels (except 
for the immediate vicinity of Idmp = 0 where any bound 
state must occur close to a gap edge and is thus hard 
to detect in a finite-size system). Hence, the criterion 
of Ref. 6 suggests that the matter-fermion sector of the 
strained Kitaev model displays non-trivial topology. 
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